A note on fractional linear pure birth and pure 
death processes in epidemic models 



Roberto Garra l5 Federico Polito 2 

(1) - Dipartimento di Fisica, "Sapienza" Universita di Roma 
t— I Piazzale Aldo Moro 5, 00185 Rome, Italy. 

Email address: rolinipame@yahoo.it 
£S) (2) - Dipartimento di Matematica, Universita degli studi di Roma "Tor Vergata" 

(3JT) Via della Ricerca Scientifica 1, 00133, Rome, Italy. 

Tel: +39-06-2020568, fax: +39-06-20434631 
Email address: polito@nestor.uniroma2.it 

ON 

(N August 30, 2011 

p I Abstract 

— In this note we highlight the role of fractional linear birth and linear death 

processes recently studied in Orsingher et al. [9] and Orsingher and Polito 
]8l, in relation to epidemic models with empirical power law distribution 

i— —i of the events. Taking inspiration from a formal analogy between the 

equation of self consistency of the epidemic type aftershock sequences 
(ETAS) model, and the fractional differential equation describing the mean 
value of fractional linear growth processes, we show some interesting 
applications of fractional modelling to study ab initio epidemic processes 

IT") without the assumption of any empirical distribution. We also show that, 

in the frame of fractional modelling, subcritical regimes can be linked 

00 t0 linear fractional death processes and supercritical regimes to linear 

fractional birth processes. 

Moreover we discuss a simple toy model to underline the possible 
application of these stochastic growth models to more general epidemic 
phenomena such as tumoral growth. 

Keywords: ETAS model, fractional branching, birth process, death process, 
5— i Mittag-Leffler functions, Wiener-Hopf integral. 



1 Introduction 

In recent years there has been a growing interest in fractional calculus modelling 
in different fields of applied sciences, from rheology to biology (see for example 
Podlubny HlOj], Debnath [g]], Mainardi [6], Diethelm [3]). It is well known that 
fractional derivatives are a good instrument to handle memory mechanisms, 
which are very useful in studying anomalous diffusion processes. Indeed, the 
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fractional derivative in the sense of Caputo is defined by a convolution of a 
power law with the ordinary derivative of the function, thus underlining its 
importance as an instrument to describe processes with emerging power law 
distributions. Recent papers []£l [Sj describe fractional processes of death and 
birth in which the integer-order derivative in the difference-differential equations 
governing the state probabilities of the related classical models, is replaced with 
the Caputo fractional derivative. In this note we apply this generalisation to 
the study of epidemic processes. In particular, we take inspiration from an 
interesting analogy between the self consistency equation of the ETAS (epidemic 
type aftershock sequences) model, used with success in statistical seismology, 
and the equation governing the mean behaviour of the fractional linear death 
model. The ETAS model, based on epidemic branching processes, llllll . is used 
to describe seismic sequences of shocks in the context of a stochastic model 
in which each aftershock may trigger other aftershocks. This model is based 
on empirical power law distributions emerging from the analysis of seismic 
catalogues, i.e., the modified Omori law on the time distribution of aftershocks 
and the Gutenberg-Richter law for the distribution of their magnitudes. These 
types of models are greatly studied in literature (see for example Sornette and 
Sornette [12] , Helmstetter and Sornette [4]) and used in assessing seismic 
risk []5]] . The self consistency equation of the ETAS model II12II is the classical 
Wiener-Hopf integral equation, but with simple handling it can be considered 
to be a fractional integral equation similar to the equation governing the mean 
behaviour of the fractional linear death model. Using this analogy we suggest 
an interpretation of the fractional death and birth processes in the framework 
of epidemic models with power law distributions (e.g. ETAS). The organization 
of this paper is the following: in Section [2] we recall the fundamental results 
on fractional death and fractional birth models; in Section [3] we describe the 
analytic formulation of the ETAS model from which we took inspiration for 
the discussion which follows in Section |4j that is the relations between the 
ETAS and fractional processes; in Section [5] in order to highlight the utility of 
this stochastic framework, we give an example of application of this formalism 
to describe tumoral growth; finally, in Section [6] we discuss the main results 
obtained. 



2 Fractional linear pure death and pure birth processes 

Recently in Orsingher et al. [ 9 ] , section 2, the fractional linear death process was 
introduced and studied. The fractional linear death process is a generalisation 
of the well-known classical linear death process. Fractionality is obtained by 
replacing the integer-order derivative in the difference-differential equations 
governing the state probabilities of the classical model with the Caputo fractional 
derivative ( |3.7| ) of order v e (0, 1). We exclude the case v = 1 as it is trivial and 
coincides with the classical case. Let n be the size of the population at time 
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and let M v (t), t > 0, v e (0, 1), be the number of components in the population 
in the fractional linear death process at time t. 

The state probabilities Pr{M v (t) = fc|M v (0) = n }, < k < n , solve the 
following Cauchy problem: 



Vp fc (t) = fj,(k + l)p fc+1 (t) - fxkp k (tl < fc < n , v e (0, 1), 

1 iffc = n , (2.1) 
if < n , 



where pi > is the death intensity and in which we consider that p„ +i(t) = 0. 
We recall the following fundamental theorem []9j page 73]: 

Theorem 2.1. The distribution of the fractional linear death process M v (t), t > 
with n initial individuals and death rates /JL k = /ik, is given by: 

pl(t) = Pr{M v (t) = k\M v (0) = n } (2.2) 

where < k < n , t > 0, v e (0, 1). The function £ v>1 (x) is the Mittag-Leffler 
function defined as 

+oo x k 

In the same paper it is also shown that the mean value EM v (t), t > 0, of the 
fractional linear pure death process satisfies the fractional differential equation 

(D v EM v (t) = -/xEM v (t), t>0,ve(0,l), 
\EM v (0) = n , 

which can be easily solved by means of the Laplace transform, obtaining 

EM r (t) = n E v l (-nt v ), t > 0, V S (0, 1). (2.5) 

In Orsingher and Polito Q8Q a the fractional pure birth process 9t v (t), t > 0, 
v e (0, 1), was studied with similar methods to those used for the fractional 
death. For the sake of our discussion, we recall only the fractional differential 
equation which governs the mean value E9T (t), t > 0, referring to the original 
paper for further details: 

fDW(t) = rEW'(t), t>0, r >0,ve(0,l), 
\ (2.6) 
EW(0) = 1, 
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where y is the birth intensity. The solution to Q2.6D reads: 



E9f(t) = £ V;1 ( r t v ), t > 0. (2.7) 

Recalling the definition of Mittag-Leffler function, we remark that asymptotically, 
it has a power law decay (see for example Ogata U71D. 

Finally, a subordination relation is verified for the fractional linear death and 
fractional linear birth processes. For example, for the fractional death, we have 
that 

M v (t) = M(v;'), t>0,ve(0,l), (2.8) 

where the equality is intended for the one-dimensional distribution and V t v , 
t > 0, v e (0, 1), is the right-inverse process of the v-stable subordinator S"L, 
x > 0, v e (0, 1), i.e. 

V t v = inf{x > 0: S* > t}. (2.9) 

3 Analytical formulation of the ETAS model 

In this section we follow Sornette and Sornette [12] to define the master 
self consistency equation of the ETAS model. We do not discuss in detail the 
modellistic aspects of the ETAS model (see for example Helmstetter and Sornette 
H4ll). Briefly, it can be considered as a simple branching model in which each 
afteshock may trigger other succeding aftershocks in a cascade process. It is 
therefore natural to regard this as an epidemic process in which a given parent 
event of a certain magnitude, occurring at time t = 0, gives birth to other child 
events with rate: 

0(t)dt = (l-fc)6>tS^H(t-r o ), t>0, (3.1) 

where H(x), x e R, is the Heaviside step function, and (1 — fc)tg/t 1+e is the 
modified local Omori law on the time occurrence of the aftershocks. In the 
following, we focus our attention only on time behaviour and decay of this 
cascade process. The self consistency equation, describing the rate of seismicity 
at a given time t, reads II12II : 

rt-t 



N(t) = (i-fc) 

Jo 

■+0O 



jV(T)0^(t-T)- (1+9) dT, (3.2) 



where J Q °° <^>(t)dt = (1 — k), is the branching ratio, i.e., the average number 
of aftershocks generated by each event; t is a time delay with respect to the 
mainshock that occurs at time zero. Note that equation Q3.2P is the classical 
Wiener-Hopf integral equation with a power law kernel. In Sornette and 
Sornette [12], and Helmstetter and Sornette [01, an interesting discussion of 
the analytic solution to this equation was given and, in particular, three different 
regimes for the seismic rate iV(t), were found: 
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• k > and 9 > 0: subcritical regime, less than one child per parent; 



k < and 9 > 0: supercritical regime. A transition from an Omori decay 
with exponent p = 1 — 9 to an explosive increase of the seismicity rate; 

9 < and k > 0: a transition from an Omori law with exponent 1 — \ 9\ to 
an exponentially increasing seismicity rate for large values of time. 



In Sornette and Sornette ||12||. a first discussion about the analytic solution of 
the self consistency equation was given. The authors found a characteristic time 
t*, function of all the ETAS parameters, acting as a cut-off between the small 
and the large time behaviour. They discussed this important feature both for the 
subcritical and supercitical regimes. For the subcritical regime, this means the 
transition from an Omori exponent p = 1 — 9 for t < t* to p = 1 + 9 for t > t* , 
i.e. it implies that the seismic rate decays slowly for small values of t. For the 
supercritical regime the authors found a transition from an Omori power law 
decay with exponent p = 1 — 9 for t < t* to an explosive exponential increase. 
In a successive paper, Helmstetter and Sornette [4] generalised their analysis, 
by taking into consideration also the Gutenberg-Richter distribution. What is 
more interesting for the sake of our discussion, is that they also analysed the 
third regime, i.e. 9 < 0. This case needs more attention, because the integral 
J ^>(t)dt ~ j Q l/(t + l) 1+9 dt becomes unbounded, thus implying an infinite 
branching ratio. The authors explained this, apparently meaningless result say- 
ing that the number of children created beyond any time t exceeds the numbers 
of children created until time t. This interpretation is also confirmed by the 
discussion on the analytic solution related to this third regime, in which another 
charateristic time t of transition from an anomalous slow decay for t < t and 
an explosive exponential increase for t > t, similar to the supercritical case, 
was found. Therefore, this case, and the supercritical case, have a modellistic 
utility proven to be useful to understand phenomena with anomalous slow 
seismic decay for small values of time, although they appear meaningless for 
large values of time. We refer to Helmstetter and Sornette [4] for a thorough 
discussion on the meaning of the regimes of the ETAS model. 
Before going ahead, we give some mathematical remarks on equation Q3.1J ). For 
example, in the last case (9 < 0), we can rewrite it in the following manner: 



N(r) = —(1 — fc)|0| 

Jo 



AT(T)t*?(r - T )- (1H9|) dT, t > 0, (3.3) 



where we have done the approximation t — t ~ t, the time delay being small 
with respect to the total time interval. By recalling the definition of Riemann- 
Liouville fractional integral of order a > 1 1 1 Oil 
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J a f(t) 

r(a) 



/(T)(t-T)- (1 " a) dT, a>0, t>0, (3.4) 
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with T(a) the Euler Gamma function, it is easy to rewrite the equation Q3.3J ) as 
a fractional integral equation. Indeed we have that 

N(t) = -M l6l N(t), t>0,9<0, (3.5) 

where A = (1 — fc)|0|r(|0|). Finally, we arrive at the following fractional 
differential equation: 

D w N(t) = -XN(t), t>0, 9<0, (3.6) 
where D' e ' is the Caputo fractional derivative of order |0| > 0, defined as 1 1 1 Oil 



p t CL 



D |e| iV(t) 



r(m-|0|) 



Mr) 



dT "' dr, (3.7) 



-m+l 



where m = \\6\~\. 

Equation ( |3.5| ) can thus be written as a simple linear fractional differential 
equation. Equation ( |3.6| ), in the case 9 e (—1,0), with the initial condition 
IV(0) = n , n e N (meaning that the initial rate of seismicity is n ), can be 
easily solved as we have shown in the previous section. 

Finally, we realize that the ETAS self consistency equation is formally similar 
to equation fl2.4| ) governing the mean behaviour of the fractional linear death 
process. Furthermore, the supercritical regime with explosive growth can be 
linked in the same way to the fractional linear pure birth process. In the next 
section, we discuss similarities and differences between these models. 



4 Relations between self consistency equation of the 
ETAS model and fractional birth or death processes 

Taking inspiration from the formal analogy shown in Section [3] here, we discuss 
the analytical behaviour of the rate iV(t), t > 0, in the framework of fractional 
death and fractional birth processes, considering, on the one hand, the behaviour 
predicted by them, and, on the other hand, the physical meaning of these 
predictions in the framework of the ETAS model. In short, the rate of seismicity 
iV(t), t > 0, is considered here as the mean behaviour of an underlying fractional 
linear death or birth process. 

Similarly to the ETAS model, starting from equation ( |3.6| ) describing the rate of 
an epidemic process, we can distinguish two regimes: 

• if X > 0, we have 

D v jV(t) = -AJV(t), t > 0, v e (0, 1). (4.1) 

As we have shown, it is identical to the fractional differential equation 
arising in the fractional death process. Physically it corresponds to the 
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Figure 1: The mean value of the linear fractional death process with v = 1 (red), 
v = 0.7 (black) and v = 0.4 (green). 

subcritical regime in the framework of the ETAS model. However we must 
underline some differences, remembering that our formal analogy is not 
rigorously an identity. We observe a strong exponential decay for small 
values of t and an asymptotic power law decay when t — > oo. Thus, we 
have a reasonable physical picture for cases with a slow decay for large 
times, but this prediction is slightly different from that of ETAS. 

In Fig. [TJ the solution to the Cauchy problem composed by ( |4.1| ) and 
the initial condition N(0) = n = 1, for some values of v, is shown; it 
corresponds to the mean value of the linear fractional death process. Note 
the fast decay for small values of t. 

• if A < we have the following equation: 

D v N(t) = \X\N(t), ve(0,l). (4.2) 

This corresponds to the fractional differential equation arising in the 
linear fractional birth process. This case is equivalent to the supercritical 
regime of the ETAS model. Indeed, we have an explosive growth. Fig. [2] 
presents the solution to the Cauchy problem composed by ( |4.2D and the 
initial condition N(0) = 1. We have an exponential growth, as for the 
supercritical regime of the ETAS model. 
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Figure 2: The mean value of the linear fractional birth process with v = 1 (red), 
v = 0.7 (black) and v = 0.4 (green). 



5 An application of the fractional pure birth process: 
theoretical model of tumoral growth 

Beyond the formal analogy with the ETAS model we can explain the utility of 
our proposal with a simple toy model about tumoral growth. In the following 
we show that the fractional linear pure birth process can give an interesting 
framework to modelling this kind of phenomena and have good agreement with 
some results in this field. Assume that we have an initial number of infected cells 
W(0) = n . We are interested in building a theoretical model about the growth 
of cancer cells with respect to time in the dynamics of metastasis. Suppose that 
this process is described by a linear fractional pure birth model. Therefore we 
have (see Orsingher and Polito [8]) that the distribution of the population size 
of cancer cells at time t is given by: 

Pm v (t) = k + n \m v (0) = n } (5.1) 
= ("° + ^ ~ ^ 2 (r) (-D%,i(-("o + r) r t r ), t > 0, v € (0, 1], y > 0. 



Now, from fl5.1[ ), we can infer the probability of a new offspring at the beginning 



of the process: 



Pr{0I v (dt) = l + n |91 v (0) = n } (5.2) 
= n [E v>1 (-n r(dtr)-E V)a (-r(n + l)(dtr)] ~"o r( ^ + 1) , 
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by writing only the lower order terms. This result shows that the probability 
of a new offspring is proportional to the time interval (dt) v and to the initial 
number of progenitors. This is an interesting picture that can be realistic for a 
large number of complex growth dynamics with emerging power law behaviour. 
From the previous analysis, we also have that the average expansion of the 
population of cancer cells is given by: 

EOT(t) = n E v>1 ( r t v ), t > 0. (5.3) 

In conclusion we have a theoretical framework that can be used to describe 
growth processes exhibiting power law expansions. 

We discussed this model by referring to the cancer growth dynamics because 
several works (see for example Dattoli et al. [HJ, West et al. II13II ") studied a 
similar dynamic starting from the Kleiber law. It is not the purpose of this 
paper going inside this vaste field of study, but we remark that we have a 
conceptual well posed model without empirical assumptions. Thus we obtain a 
probabilistic point of view that is more fundamental with respect to the dynamic 
of these processes. Finally, we must observe that, in order to make the model 
more realistic, it is possible to directly generalise the fractional birth process by 
introducing a saturation threshold. This will be matter of a future paper. 
We also notice that the role played by the allometric coefficients in the Kleiber 
law, in our model, is played by the real order of derivation v. These allometric 
coefficients are empirically established; we have the same picture but with only 
one free fitting parameter. Beyond the probabilistic clear view of the processes, 
this is also a strong point of this theory, 

6 Discussion 

We now move to understand the differencess between the ETAS model and 
the fractional epidemic processes predictions. First, we must identify the real 
order of derivation v with the parameter 6 of the Omori law and A with 
X - (1 - fc)|0|r(|0|) of the ETAS model. Although we do not have exactly 
the same cut-off between different regimes as in the ETAS model, we obtain 
some reasonable agreement for various ranges of the parameters. Moreover, in 
Helmstetter and Sornette PH , by discussing the self-consistency equation for the 
supercritical regime, a cut-off between slow decay for small times and explosive 
growth was found. In our case, the supercritical regime is less meaningful, 
presenting immediately an exponential growth. On the other hand, from the 
self-consistency equation fl3.2J ) and by considering the physical behaviour in the 
subcritical regime, the agreement of the two models is clear, both describing 
epidemic processes. 

At this point, by keeping in mind the relations between the ETAS model and the 
fractional linear birth or death processes, we make some remarks. We have a 
generalised analytic model with clear relations with epidemic models. We also 
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notice that in many epidemic processes in complex systems (as in seismology 
and biology), empirical power law behaviour emerges. We have seen that the 
prediction on the rate of fractional processes, i.e., the Mittag-Leffler function 
(generalised exponential), takes into account this behaviour asymptotically. We 
conclude that the fractional processes we treated in this paper, are an important 
instrument in order to study epidemic processes which exhibit a rapid decay for 
small times and slow asymptotic power law decay for large times. 
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